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We calculate the yields of quarkonia in heavy ion collisions at RHIC and the LHC as a function 
of the transverse momentum. Based upon heavy quark effective theory, our results include both 

$H ' color-singlet and color-octet contributions and feed-down effects from excited states. In reactions 

with ultra-relativistic nuclei, we focus on the consistent implementation of dynamically calculated 
nuclear matter effects, such as coherent power corrections, cold nuclear matter energy loss and 
the Cronin effect, in the initial state and collisional dissociation of quarkonia in the final state, as 
they traverse through the QGP. Theoretical results are presented for J/ip and T and compared to 
experimental data where applicable. At RHIC, a good description of the high-pr J/ip modification 
observed in central Cu+Cu and Au+Au collisions can be achieved within the model uncertainties. 
We find that J/ip measurements in proton-nucleus reactions are needed to constrain the magnitude 
of cold nuclear matter effects. At the LHC, a good description of the experimental data can be 

CIh 1 achieved only in mid-central and peripheral Pb+Pb collisions. The large five-fold suppression of 

prompt J/ip in the most central nuclear reactions may indicate for the first time possible thermal 
i-C J effects at the level of the quarkonium wavefunction at large transverse momentum. 
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I. INTRODUCTION 

Melting of heavy quarkonium states, like the J/t/j and the T, due to color screening in a deconfmed quark-gluon 
plasma (QGP) (l| has been proposed as one of the principal signatures for its formation. An expected experimental 
consequence of this melting in the thermal medium created in heavy ion collisions (HIC) is a suppression of the yield 
of heavy mesons, when compared to their yield in nucleon-nucleon (NN) collisions scaled with the number of binary 
interactions. 

In a simplified picture, one can think of a quarkonium as a QQ pair in a color-singlet bound state, where the heavy 
quark (Q) and the anti-quark (Q) are separated by distances ~ l/{m,Qv) smaller than 1/Aqcd- Here, v ~ a s {rriQv) 
is the relative velocity between Q and Q, and the state is held together by an effective potential interaction. The sizes 
of the higher excited states are larger, and the sizes of the bottomonia are smaller than the corresponding charmonia. 
The presence of a thermal medium screens the interactions between the Q and Q and leads to melting at some 
characteristic temperature [2| that depends on the meson. This picture suggests that it may be possible to observe 
sequential melting of narrower quarkonia as we explore thermal media of increasing temperatures 3]. Several studies 
(see [4j for a recent review) have calculated the modification of quarkonium yields in collisions at SPS, RHIC, and 
theLBC grllil. 

A more sophisticated description of quarkonia is provided by Heavy Quark Effective Theory (HQET) [Tl( . In this 
picture, the QQ state in the color-singlet combination is the lowest order Fock space component of the full quarkonium 
wavefunction. Higher Fock components in the wavefunction have additional partons and are suppressed by positive 
powers of the small parameter v. In our calculation we will take into account the Fock components that allow for both 
color-singlet and color-octet QQ contributions. The HQET formalism has been used to calculate differential yields of 
heavy mesons as a function of the transverse momentum, pt, in p+p collisions [Ill4l5|. The picture of the production 
process in HQET is as follows. The initial hard collision produces a "proto-quarkonium" , a short distance (~ 1/ttiq) 
QQ pair where the Q and the Q can be in the color-singlet or the octet state. Since this process is short-distance, 
it can be calculated using perturbative QCD [lf|. This "proto-quarkonium" then evolves into the quarkonium state 
with probabilities that are given by non-perturbative matrix elements. For color-octet states, this evolution process 
also involves the emission of soft partons to form a net color-singlet object. The formed final state has a hierarchy of 
Fock space components that are given by HQET scaling rules. In this paper, we combine the HQET formalism with 
cold nuclear matter (CNM) effects and the effects of propagation of quarkonia through a hot QGP to calculate the 
final yields in HIC. We then obtain the nuclear modification ratio, 

der^/dydpr 

HabKPt) ; q3~ , (1) 

N AB d(J PP '/dydPT 

in reactions with heavy nuclei. We present results for Au+Au and Cu+Cu collisions at RHIC at \/~S = 0.2 TeV per 
nucleon pair and Pb+Pb collisions at LHC at s/S = 2.76 TeV per nucleon pair. 

We emphasize upfront that our goal is not to fit data in heavy ion reactions but to investigate systematically the 
relative significance CNM and QGP effects and discuss ways in which these effects can be better constrained. Our 
formalism is limited to high transverse momentum and not applicable to pT-integrated yield. The CNM [17| effects 
we consider include nuclear shadowing (here implemented as coherent power corrections) and initial state energy loss. 
These are calculated in the same framework previously used for light hadrons [l8l - [24l ] and open heavy flavor [2a, [26j ■ 
We also consider transverse momentum broadening [18|, [2|| [29| (as a model for the Cronin effect). We give results for 
the p(d)+A collisions at RHIC and the LHC with and without Cronin. Since the CNM effects are generated on a short 
time scale ~ (lfm/c)^4 1 / 3 /7 (the Lorentz contracted size of the nucleus in the lab frame) compared to the formation 
time of the quarkonia, we take the QQ pair in the final state in the CNM calculation as the "proto-quarkonium" 
state. 



We include the effects of QQ propagation through the thermal medium by solving the rate equation [33 . |33| 
describing the change in the production yield of quarkonia as a function of time. The microscopic input for the 
equations arc the time scales for the formation and the dissociation of the quarkonia. Since the size of proto- 
quarkonium is small, it is typically assumed that the formation time if or m for the "proto-quarkonium" to "expand" 
into the quarkonium wavefunction, is not significantly modified by the presence of the QGP [8(. Hence, t{ orm in the 
rest frame of the meson is ~ \ / (rriQV 2 ) . We estimate tf 01 m by dividing the radial size of the wavefunction with the 
radial velocity and assign a factor of 2 uncertainty to it. (Some equilibrium approaches [5| include effects of a finite 
time of formation by suppressing the decay rate.) 

For dissociation times, most studies use results for thermally equilibrated quarkonia [a, 0, Q- For a quarkonium 
stationary or slowly moving through the thermal medium it might be safe to assume that it reaches thermal equilibrium 
with the medium. In this picture, the wavefunction of the quarkonium (the lowest Fock component) is better described 
by the solution of the Schrodinger equation with a modified potential that depends on the temperature [2|, M-|S | . The 
potentials used have a firm footing based on lattice calculations of the potential between two static charges [41 1. The 
imaginary part of the potentials have been calculated in [30,[3l|. Other approaches use a T— matrix formalism [7| to 
calculate decay rates of thermal states. 

A narrow object produced at t ~ and moving rapidly through the QGP may not have sufficient time to thermally 
equilibrate with the medium. As it propagates through the hot and dense matter, collisions with thermal gluons 
can dissociate the quarkonium on a time scale tdiss [32h35| . In this paper, we take this high-pr limit and explore 
the consequence of assuming that the wavefunction of the meson is the same as in the vacuum, and that the main 
impact of the medium is to dissociate the meson as it propagates through it. This dissociation mechanism has 
been phcnomcnologically successful at describing the modification of the yields of open heavy flavor mesons as a 
function of the transverse momentum pt- At RHIC, this mechanism is consistent with the large suppression of non- 
photonic electrons I37H39J . The predictions [33j for D-meson suppression at the LHC arc compatible with preliminary 
experimental results [401 ] . The formation time for quarkonia is also small, like it is for the heavy- light mesons. 
Therefore, in-medium formation and dissociation can modify the quarkonia yields. We find the dissociation time idiss 
for the color-singlet Fock component of the formed quarkonium wavefunction by calculating the rate of broadening 
of the wavefunction due to in- medium interactions, as described for open heavy flavor in [32J. The color-octet state 
must emit soft gluons in order to map onto the color neutral quarkonium. Since soft gluon emission occurs on a short 
time scale oc lu and does not change the QQ kinematics, the corresponding Fock component can be thought of as two 
"effective" heavy quarks, consisting of the ones produced in the hard collision and the emitted soft partons. In this 
approximation, the broadening of the wavefunction component onto which the short distance color-octet QQ map is 
proportional to the broadening of color-singlet state and gives similar dissociation rates. One important difference 
with the heavy-light mesons is that the probability of reformation after dissociation is small since the probability 
of refragmentation to quarkonia is small. In principle, recombination [42( with thermal Q and Q from the medium 
can re-generate quarkonia. Since the population of thermal heavy quarks is exponentially suppressed (~ e~ mT ' T ), 
especially at large pr, we ignore this process. 

Our results for the Raa of J /if) mesons, obtained using this formalism, are consistent with PHENIX and STAR data 
for central Au+Au and Cu+Cu collisions. p+A quarkonium measurements at high pt are necessary to accurately 
constrain the magnitude of CNM effects. Results from ATLAS and CMS at the LHC for the J/ip suppression are well 
reproduced by the calculation only in peripheral and mid-central collisions. It is difficult (with or without transverse 
momentum broadening) to obtain a five- fold suppression of prompt J/ip, as reported by the CMS experiment at 
the LHC, without assuming thermal effects at the level of the quarkonium wavefunction. For the Ts, this approach 
predicts that the suppression for the T(2S*) and T(35*) at high px is comparable to the T(15) suppression. This 
prediction is also different from the expectation of the equilibrium models, which predict higher suppression for the 
excited states. Differential high transverse momentum data is also necessary to shed light on the possible thermal 
modification of boosted Ts at the LHC. 

Our paper is organized as follows. In Scction[H]we present the calculation of the p+p baseline yields of quarkonia. In 
Section Hm we discuss how CNM effects modify the production yields of quarkonia in nuclear collisions. In Section HVl 
we evaluate the wavefunctions needed for the evaluation of the formation and dissociation timescales. In Section [V] 
we calculate the dissociation time scale in the QGP and discuss the dynamics of quarkonia as they propagate through 
the strongly-interacting medium. Results for the nuclear modification factor in p+A and A+A collisions are given 
in Section IVII Wc present our conclusions in Section IVIII Wc illustrate the details that go into the evaluation of the 
baseline quarkonium production in Appendix [X] 



II. QUARKONIUM PRODUCTION IN P+P COLLISIONS 

In this section we describe the production of quarkonia at high transverse momentum in p+p collisions. It provides 
a baseline for the calculation of the nuclear modification factor Rab(pt) defined above in Eq. [T] It also gives the 
initial unquenched spectrum of "proto-quarkonia" in the heavy ion collision, excluding CNM effects and effects of the 
propagation of the quarkonium states through the QGP medium. 

The dominant processes in evaluating the differential yield of heavy mesons as a function of px the are 2 — > 2 
processes of the kind g + q — >• H + q, q + q — » H + g and g + g — >• H + g, where H refers to the heavy meson. We label 
the process generically as a + b — > c + d, where a and b are light incident partons, c refers to H and d is a light final 
parton. Given the matrix element for the process, M a b^cd, the cross section has the form 

dx a cj) a {x a , n f )(/) b (x b , Hf) — m T v XgX b -j(ab ->■ cd) , (2) 



dprdy J — — ^-— ^ Xa _m| e , - - d{ 



where momentum-energy conservation fix 



1 x a VSm T e y - m 2 H 
VS XaVS — rriTe y 



and we take the factorization and renormalization scales /x/ , /i^ to be rriT = \Jv\- + m # • The invariant cross section 
is given by 

da = \Ml_ 

dt 167TS 2 - l ' 

A. HQET calculation of da/di 

We use Heavy Quark Effective Theory (HQET) [llj results to calculate the production of quarkonia in p+p collisions. 
HQET provides a systematic procedure to compute any quantity as an expansion in the relative velocity v of the 
heavy quarks in the meson. For example, the wavefunction of the J/ip meson (analogous expressions hold for the 
ip{2S), T, T(25) and T(35)) is written as 

I J/i>) HQQ([ 3 Si]i)> + Oiv^QQ^Sohg)) + Oiv^lQQtfSJsgg)) 

+ Oiv^lQQtfPajsg)) + O^QQ^P^g)) + 0(v 1 )\QQ([ 3 P 2 ] s g)) + ■■■ 

The differential cross section for the prompt (as opposed to inclusive which includes contributions from B— hadron 
decay) and direct (as opposed to the indirect from the decay of heavier charmed mesons) production of J/tp can also 
be calculated in HQET. It can be written as the sum of the contributions, 

da(J/i>) = da(QQ([ 3 5 1 ] 1 ))(0(QQ([ 3 5 1 ] 1 ) -► J/VO) +da(QQ([ 1 So}s))(0(QQ([ 1 So]s) -+ J/V)) 

+ rfa(QQ([ 3 5 1 ] 8 ))(C(QQ([ 3 5 1 ] 8 ) -> J/j,)) + da(QQ([ 3 P } s ))(O(QQ([ 3 P Q h) -> J/V)) (6) 

+ rfa(QQ([ 3 P 1 ] 8 ))(G(QQ([ 3 P 1 ] 8 ) -> J/J,)) + da(QQ([ 3 P 2 ] s ))(0(QQ([ 3 P 2 ] s ) -> J/j,)) + ■■■, 

where the quantity in the brackets [ ] represents the angular momentum quantum numbers of the QQ pair in the Fock 
expansion. The subscript on [ ] refers to the color structure of the QQ pair, 1 being color-singlet and 8 being color- 
octet. The dots represent terms which contribute at higher powers of v. The short distance cross sections da(QQ) 
correspond to the production of a QQ pair in a particular color and spin configuration, while the long distance matrix 
clement (0(QQ) — > J/ip) corresponds to the probability of the QQ state to convert to the quarkonium wavefunction. 
This probability includes any necessary prompt emission of soft gluons to prepare a color neutral system that matches 
onto the corresponding Fock component of the quarkonium wavefunction. 

Power counting rules tell us that contributions from the color-octet matrix elements in Eq. ([6]) are suppressed by 
v 4 compared to the color singlet matrix elements. More specifically, 

(0(QQ([ 3 5i]i) -► J/VO) = 0(m 3 Q v 3 ) , 
(O(Q0([ 3 ^] 8 ) -4 J/VO) - 0{m%v 7 ) , 
(OiQQ^Soh) -> J/VO) - 0{m%v 7 ) , 
(0(QQ([ 3 Pj] s ) -> J/iP)) = 0{m%v 7 ) . 



These operators are multiplied by the short distance differential cross sections, which are related to the probability 
to create QQ pairs in specific quantum states. Since these are short distance operators, these can be calculated in 
perturbation theory. We use the expressions for the short distance color-singlet cross sections given in [4J, |44j and 
the color-octet cross sections given in jl^ - 4l4j ] . 

The case of the P-wave bound states (xcOi Xci and \c2, sometimes collectively referred to as XcJ and the corre- 
sponding states of the b quark) is slightly different. The wavefunction of Xc states can be written as, 

ha) HQQ([ 3 P,/]i)> + 0(«)|QQ([ 3 Si] 8 <?)) + O(v 2 )\QQ([ 1 S ]sg)) + 0(v)\QQ([ 3 Dj] 8 g)) 
+ Oiv^lQQ^P^sg)) + 0(v 2 )\QQ([ 3 Pj] 8 gg)) + ■■■ 

The color-singlet state QQ[ 3 Pj]i and the color-octet state QQ[ 3 Si]$ contribute to the same order in v because of the 
angular momentum barrier, and hence both need to be included for a consistent calculation in v. For the calculation 
of the production cross section, we consistently take the contributions to the lowest order in v as follows. Xc as, 

da( Xc j) = do(QQ{[ 3 P J ] 1 )){0{QQ([ 3 Pj] l ) -> XcJ )} + da(QQ([ 3 Si] 8 ))(0(QQ([ 3 Si] s ) -> XcJ )) + • • • (9) 

Similar expressions hold for the xt,(lP), Xb{2P) and Xfc(3P) mesons. The expressions for the short distance coefficients 
are given in [131 ] . The scaling of the matrix elements are given as follows, 

(0(QQ([ 3 P J ] 1 ) -► xci)) = 0(m 5 Q v 5 ) , 

(0(QQ([ 3 5i] 8 ) -> Xci)) = 0(m 5 Q v 5 ) . 

We now look at the values of the color-singlet and color-octet expectation values and use them to give theoretical 
results for the production of quarkonia. 

B. Matrix elements for QQ production 

The scaling of the operators is described in Eqs. (J7J, (fTO)) but their exact values can not be calculated using 
perturbation theory. In this work, following [12, LL3| we use the values of color-singlet operators calculated using the 
potential model. The expressions and the values for the color-singlet operators are given in [12|, [lj, |45| . The values 
we obtain by solving the non-relativistic wavcf unctions: 

(0(cc([ 3 5i]i) -> J(n = 1))) = 3iV c |fl - 2 f )|2 = 1.2(GeV 3 ) , 

(O(cetfS ]i) -* J(n = 1))) = | x (0(cc([ 3 5i]i) -► J(n = 1))) 

(0(cc([ 3 Po]i) -> Xco(n = 1))> =3JV e |fl "=; (0)|2 = 1.63 x lO-^^GeV 3 ) , ( u ) 

(0(cc([ 3 Pi]i) -► Xco(n = 1))> = 3 x (0(cc([ 3 Po]i) -► X co(" = 1))) , 

(C(cc([ 3 P 2 ] 1 ) -> X co(n = 1))> = 5 x (O(cc([ 3 P ]i) -> X co(n = 1))) , 

where R(0) is the radial wavefunction at the origin, R'(0) is the first derivative of the radial wavefunction at the 
origin, and n refers to the radial quantum number. A slightly better fit to the data is given by the following values 
of the matrix elements (within ~ 20% of the ones quoted above), 

(0(cc([ 3 5i]i) -> J(n = 1))) = 1.5(GeV 3 ) , 

(0(cc([ 3 P a } 1 )^Xco(n = l))) =2.04xlO- 1 m2( Ge V 3 ). [LZ > 

and we use these for the final differential cross sections. 

The values of the color-singlet operators for bottomonia are given in [14]], which we reproduce here: 

= 10.9(GcV 3 ) , 
= 1.01 x 10- x mg(GeV 3 ) , 

= 4.5(GeV 3 ) , (13) 



(0(bb([ 3 S 1 ] 1 )^T(n=l))) =3Nj-^M 

(0(bb([ 3 Po]i)^Xw(n = l))) = 3A C ^%^ 

(0(bb([ 3 S 1 ] 1 )^T(n = 2))) =3N c l -^2l 

(O(bb([ 3 P ] 1 )^ Xb0 (n = 2))) = 3N C &^ 

(0(bb([ 3 S 1 } 1 )^T(n = 3))) =3N C ^^1 
A better fit to data is given by: 



3.64 x 10- 2 m2(GeV 3 ) , 
= 4.3(GeV 3 ) . 



(C(66([ 3 5i]i) -► T(n = 1))) = 11.99(GeV 3 ) , 

(O(66([ 3 P ]i) -► Xbo(n = 1))) = 1.01 X lO^mKGeV 3 ) , 

(C(66([ 3 5i]i) -> T(n = 2))) = 5.40(GeV 3 ) , (14) 

(O(66([ 3 P ]i) -»■ Xbo(« = 2))) = 3.64 x 10- 2 m 2 (GcV 3 ) , 

(C(66([ 3 5i]i) -► T(n = 3))) = 5.16(GeV 3 ) . 



These values are again within 10%-20% of the originally suggested ones. 

The color-octet operators can not be related to the non-relativistic wavefunctions of QQ. Following |12| - [l4| . we fit 
them to the data. We give these values below: 



(0(cc([ 8 Si] 8 ) 
(0(cctfSo]a) 

(O(cc([ 3 P ]s) 

(0(cc(i 3 Pi} 8 ) 

(0(cc([ 3 P 2 } 8 ) 

(0(cc([ 3 S 1 ] s ) -i 

(0(cc([ 3 ^] 8 ) 

r3 Si] 8 ) 



+ J{n = 1) 

+ J{n = 1) 
+ J{n = 1) 
■+ J(n = 1) 
+ J{n = 1) 
Xco(n = 1) 



-> Xd{n = 1)) 
(0(cc([ a ^] 8 ) -+ X c2(n = 1)) 

The octet matrix elements for the bottomonia arc taken from a fit to TeVatron data and arc as follows: 



= 3.31 x l(T 3 (GeV 3 ) , 

= 7.14 x l(T 2 (GcV 3 ) , 

= -4.21 x 10" 3 m 2 (GcV 3 ) , 

= 3 x (O(cc([ 3 P ] 8 ) -> J(n = 1))) , 

= 5 x (0(cc([ 3 P o ] 8 ) ^ J(n = l))>, 

= 1.02 x l(T 3 (GeV 3 ) , 

= 3 x (0(cc([ 3 51] 8 ) ^ X co(« = 1))), 

= 5x(0(cc([ 3 51] 8 )^ X co(« = l))). 



(15) 



(0(66([ 3 5i 

(0(bb([3P Q 
(0(66([ 3 5!] 8 
(0(66([ 3 5i 
(O(6&([^o 
(O(65([3P 
(0(M([ 3 Si] 8 
(0(&6([ 3 Si 

(O(&6([3P 







>T(n = 1) 
>T(n = 1) 

> T(n = 1) 

Xbo("- = 1) 

> T(n = 2) 
>T(n = 2) 

> T(n = 2) 
) -> X6o(n = 2) 
8 )^T(n = 3) 
8 ) -> T( ? i = 3) 
8 )^T(n = 3) 



= 1.0x 10 _1 (GeV 3 ) , 
= -1.01 x 10 _1 (GeV 3 ) , 
= 1.01 x 10" 2 mg(GeV 3 ) , 

1.9 x 10" " 
= 4.18x 10" 
= -6.84x 10~ 2 (GcV 3 ) , 
= 1.37x 10" 2 mg(GeV 3 ) , 
= 1.2 x 10"" 

2.96 x 10" 



" 2 (GeV 3 ) , 
* -2 (GeV 8 ) 



(16) 



'- 2 (GeV 3 ) , 
" 2 (GcV 3 ) 



= 7.20 x 10" 2 (GeV 3 ) , 

= -1.56 x 10~ 2 to 2 (GcV 3 ) . 



Feed-down contributions 



In this project we will focus on the yields of J/ip, T(15), T(25), and T(3S) mesons. As discussed in Section Hi Bl 
we will use the production from the color-singlet and color-octet "proto-quarkonia" . This gives the direct production 
of each quarkonium species. Higher excited states of the mesons can also undergo decay to the states of lower energy. 
These contributions arc difficult to disentangle experimentally, and therefore we include the feed-down contributions 
to obtain what is called the prompt yield. For example, Xcj mesons contribute to the prompt yields of J/ip. More 
details about the feed-down contributions to J/ip are given in Appendix [AJ Finally, B hadrons can also decay to J/ip 
with a branching fraction ~ 10~ 3 . Adding their contribution gives the full inclusive yields of J/ip. In particular, at 
high pt, the contribution to the inclusive yield from the decay of B— hadrons is substantial, and can possibly even 
dominate production. At the LHC and the TeVatron the B— decay contributions have been separated, while RHIC 
reports the inclusive yields. 

In Section III Dt we show the comparison of J/ip production yields in p+p and p+p collisions to experiments and 
in Section Hi El we repeat the exercise for Ts. 



D. J/ip production in p+p and p+p collisions 

Shown in Fig.[TJare yields for JM production at the TeVatron, RHIC and the LHC at a/5 = 1.96, 0.2, and 7 TcV, 
respectively. The data from CDF [46| at the TeVatron and from ATLAS j47j at the LHC shows the prompt production 
and can, therefore, directly be compared to the theoretical curves of the prompt yields obtained from HQET and 
shown in Fig. Q] The PHENIX [48( and STAR [49| data from RHIC shows the inclusive yields. Using previous 
results for the production of B— hadrons [32[ and calculating their decay to J/ip, we have estimated the feed-down 
contribution from B— hadrons to the J/ip yields. Our model for B— decay gives results that match well with the 
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FIG. 1: (Color online) J/ip production yields multiplied by the branching ratios, B(J/ip — ¥ ee) ~ B(J/ip — > /U/u) ~ 5.93%. 
The upper left panel corresponds to data for the yield of J /ip from CDF at \/~S = 1.96 TeV [4q |. The upper right panel is 
for RHIC at yS = 0.2 TeV with inclusive data from the PHENIX experiment [481 ] (using the fit given in the paper) and the 
STAR experiment 0. The bottom panel has data from the LHC at y/S = 7 TeV from the ATLAS collaboration [U The 
darker colored error bars (dark green or red online) give the systematic errors and the lighter colored give the statistical errors 
in the experimental data. The solid line (blue online) is the theoretically calculated result. We give a 25% uncertainty band as 
a rough estimate to reflect that the color singlet and octet matrix elements are not known. The CDF results are for rapidity 
\y\ < 0.6, the RHIC and LHC results are quoted per unit rapidity at mid rapidity. ATLAS and CDF also give the feed-down 
contribution from B— decay, which we calculate and show with dashed curves (cyan online). 

observations from CDF and ATLAS. At RHIC we find that at least up to pt ~ 20 GeV, the feed-down contribution 
is smaller than the prompt production but is not insignificant. We attribute this fact partly to the hardness of the 
HQET spectrum. 

We see that the results for HQET arc a bit less steep than observed in the data. This could be due to higher 
order corrections to HQET, though NLO calculations also suffer from the same problem [50j . The results for RHIC 
at VS = 0.2 TeV also provide the baseline for the p+p yield, which we use to calculate Raa- For the LHC, we need 
the baseline at a center-of-mass energy, V5 = 2.76 TeV, which is given in the Appendix lAl 

E. T production in p+p and p+p collisions 



We now turn to the calculation of the yields of bb states. Since the lifetimes of T(15), T(25) and T(35) are all 
greater than 100s of femtometers/c (fm/c), we give the results for these separately. If only the T(1S) yields are 
measured, then we can obtain that by adding to the direct yields, the feed-down from T(25) (branching fraction 31%) 
and T(35) (branching fraction 16.4%). We see that the values of matrix elements that give a good agreement with 
the LHC data [5l[ give smaller than expected yields at the TcVatron [52[. We expect that the deviation from the 
measured yields will be smaller at yo = 2.76 TeV, for which we calculate Raa- The difference in overall normalization 
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FIG. 2: (Color online) T yields multiplied by B(T — >■ jj.fi) ~ B(T — >■ ee) ~ 2.48%. Left panel shows data from the TeVatron 
at 1.8 TeV [H2]. The right panel corresponds to data from the LHC at 7 TeV [m|. The solid line (brown) is for T(1S), the 
dashed (red) for T(2S) and the dotted (gray) for T(3S). 



cancels out in the calculation of Raa- 



III. COLD NUCLEAR MATTER EFFECTS 



In heavy ion reactions, the production yields of energetic particles are always affected by cold nuclear matter 
(CNM) effects. These include nuclear shadowing, initial state energy loss and transverse momentum broadening (also 
interpreted as the origin of the Cronin effect). In our calculation, the CNM effects are calculated from the elastic, 
inelastic and coherent scattering processes of partons in large nuclei [26[ . 

1. Nuclear shadowing: Shadowing generally refers to the suppression of the inclusive deep- inelastic scattering (DIS) 
cross section per nucleon in reactions with a nuclear target relative to the corresponding cross section in reactions 
with a proton targ et. For Bjorken x larger than 0.25 (EMC region) the effect is mainly due to collective effects 
in the nucleus 15311 . In this region we use the parametrization provided in [5_4j. For x < 0.1, power-suppressed 
resummed [20, [2l| coherent final-state scattering of the struck partons leads to suppression in the observed cross 
sections. These are included in our calculation and modify the momentum fractions of the incident partons in 
Eq. d as follows EJ: 



1 



eM 1/5 



i) 



-t- 



Xb = x b 



1 



eM 



1/3 



1) 



-U + TO;? 



(17) 



In Eqs. (|17p i, u are the relevant Mandelstam variables at the partonic level and m c , to^ are the masses of the 
struck partons. In this case we treat the octet QQ state as a massive gluon with m c = n%H- Immediately after 
the hard scattering the unexpanded color-singlet state does not couple to the medium. (Equivalently, its color 
factor is 0.) Following previous studies J33J], we use (£? A l / 3 ) q . g m (2p 2 L/\) q , g , which yields (£ 2 ) q s=s 0.12 GeV 2 
and {i 2 ) g ~ 0.27 GeV 2 . 



2. Initial state energy loss: Before the large Q 2 parton scattering process, the incoming partons lose radiatively 
a fraction of their energy due to multiple interactions in the target nucleus. If the colliding partons a, b lose 
fractional energy e a ^ = „ "' b , to satisfy the same final-state kinematics they must initially carry a larger 
fraction of the colliding hadron momentum and, correspondingly, a larger value of x. This can be implemented 
in Eq. ([2]) by the following modification 



4>a(Xa) -> 4>a(x a ) = 4>a 



4>b{Xb) -> 4>b(xb) = 4>b 



Xb 



l-e b 



Xah < 1 



(18) 



in the parton distribution functions 4> a .b/N(x a .b, [if). The medium-induced radiative corrections factorizc as a 
standard integral convolution [361 ] . If P q ^ g (e) is the probability density for quarks and gluons to lose a fraction 



e of their energy discussed above, it can be implemented in the cross section calculation as follows j24| : 

dx(f> q (x) ► [ dx [ decj) q (-^—)p q (e), I dx<j> g [x) > [ dx [ de (f> q ( -^—) P g (e) . (19) 



The calculation of P(e) is described in |22| . 

3. Cronin effect: In p+A and A+A reactions, Cronin effect can be modeled at the level of the pr-differential cross 
sections by including the fcy broadening of incoming partons that arises from initial-state scattering [18l . |29| ] . 
This involves relaxing the assumption that the incident partons a, b in Eq. ([2]) are collinear and allowing them 
to carry a transverse momentum kx a , b with model distributions f(kx a , b)- The advantage of using a simple 
normalized Gaussian form for /(k^a, b) is the additive variance property: 

(4a, b)AB = (4a, b) N N + (4a, b)lS , (4a, b)lS = (j^) £ ■ (20) 



where (k Ta D )/s accounts for increased momentum broadening in collisions with nuclei. Specifically, in Eq. (f20l 
fj? = 0.12 GeV 2 , X g = {Cf/Ca)\ = 1 fm and £ ~few is a numerical factor that accounts for the enhancement of 
the broadening coming from the power-law tails of the Moliere multiple scattering. For light final-state partons, 
the Cronin effect is well studied phenomenologically |18(. For D- meson and B- meson production the Cronin 



effect was discussed in [33j. We note that a 33% reduction of the logarithmic enhancement factor £ relative to 
the one used in the above manuscript gives a better description of the px ~few GeV light hadron production and 
this is what we use in our work. For quarkonia, the effect of initial-state transverse momentum broadening is not 
well understood. To our knowledge, this is the first attempt to calculate the Cronin effect for QQ production. 
We do not try to constrain phenomenologically its magnitude and point out that data from p+A reactions is 
needed to constrain all CNM effects. We find that including transverse momentum broadening in the same way 
as is done for light final states reduces the suppression for pr ~ 10 GeV significantly and may actually lead to 
a small enhancement of the charmonium cross section. It is not clear, given the large error bars, whether the 
Au+Au and Cu+Cu data at RHIC are better described by the Cronin calculation. In contrast, it is evident that 
the Cronin effect is not consistent with Pb+Pb data at the LHC, which sees a large attenuation at pt ~ 8 GeV. 
In Section IVl Al we will discuss the yields without the Cronin effect. We will discuss the phenomenology of the 
Cronin effect in Section IVlB I 

IV. QUARKONIUM WAVEFUNCTIONS 

The color-singlet matrix element is proportional to the square of the value of the I th derivative of the radial 
wavefunction of the QQ Fock state at the origin, where I refers to the orbital angular momentum of the bound state. 
Also, for the calculation of the collisional dissociation rate, we need the form of the wavefunctions in momentum 
space. In this section we describe the calculation of these bound state wavefunctions of the heavy mesons. In the 
calculations, we will assume that high px quarkonia do not thermalize and, therefore, we will use only the zero 
temperature wavefunctions. 

A. The instant form wavefunctions 

Making the usual separation of the radial and angular parts of the wavefunction, ip(r) = Y[ n (f)R n (r) 1 the Scrodinger 
equation for the radial part of the QQ state can be written as 



1 d 2 1(1 + 1) 
2/i dr 2 2 fir 2 



V(r) 



rR nl (r) = (E nl - 2m Q )rR nl (r) , (21) 



where p = —£■ is the reduced mass (not to be confused with the renormalization scale), I is angular quantum number, 
and n the radial quantum number, (n — 1 is the number of nodes in the position space wavefunction.) V(r) is the 
potential between the two heavy quarks, which can be estimated from the lattice. We use the potential described 
in [2| . The results are compactly written in Table U and Table [Til where we approximate the radial wavefunction as 
R(r) = N r r l Tli=i tn -x(l — -£-) exp(~ r ^ 2a >\ where \ l r are the nodes of the wavefunction in radial coordinates. 
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For the calculation of the dissociation rate, we also require the form of the wavefunction in momentum space, 

^(k) = / ' d 3 r^(r) = YT{k)R(k) ■ (22) 



ip(k) satisfies the normalization condition J A J? 3 \ip(k)\ 2 = (2 K 3 J dkk 2 \R(k)\ 2 = 1. The momentum space wave- 
functions are readily related to the light cone wavefunction as described in [23]. In order to simplify the cal- 
culations of the dissociation rates, we approximate the momentum space wavefunctions by the form R„i(k) = 
Nkk l Hi=i n -i(l — -rr)exp( -fc ^ 2k ^, where A^ are the location of the nodes in the momentum space wavefunction. 

Our results are again presented in Table [fl and Table [TTJ The color-singlet operators are given by the expressions 
Eq. CLE]). 

I n EbjGeV) a(GeV)- 1 fc 2 (GeV) 2 iV r (GeV) 3/2 iV fc (GeV)- 3/2 A^GeV)" 1 A fc (GeV) Meson 

XcO, Xcl, Xc2 

TABLE I: Charmonia wavefunctions. I refers to the angular momentum of the QQ state, while n is the radial quantum number. 
For charmonia, only n — 1 wavefunctions are important. 

I n Eb(GeV) a(GeV)" 1 fc 2 (GeV) 2 iV r (GeV) 3/2 iV fe (GeV)- 3/2 A^GeV)" 1 A fc (GeV) Meson 



1. 


0.7003 


1.829 


0.2988 


0.6071 


58.55 


1. 1. 


0.2678 


2.216 


0.2036 


0.1677 


141.3 



1. 


1.122 


1.007 


0.9854 


1.486 


2. 


0.5783 


1.446 


0.4784 


1.479 


3. 


0.2139 


1.768 


0.3199 


1.813 


1. 1. 


0.7102 


1.309 


0.5832 


0.6251 


1. 2. 


0.325 


1.588 


0.3967 


1.004 


1. 3. 


0.05109 


2.14 


0.2183 


0.5635 


TABLE II: Bottomonia wavefuncti 


ons. 1 refers to the 


number. 











23.92 






rib, T(n = 1) 


103.7 


1.304 


1.012 


rib, T(n = 2) 


236.4 


1.106, 2.751 


0.615, 1.457 


77b, T(n = 3) 


37.9 








Xbo,i,a(n = 1) 


196.9 


0, 2.122 


0, 1.105 


Xbo,i,2(n = 2) 



758.4 0, 1.814, 3.621 0, 0.694, 1.453 Xbo,i,2( n = 3 ) 

I refers to the angular momentum of the QQ state, while n is the radial quantum 

Note that, as mentioned earlier, the QQ pair produced in the hard collision in a color-octet state has to emit 
soft gluons to overlap with the wavefunction of the color neutral quarkonium. These soft gluons do not change the 
kinematics of the QQ but simply ensure that the final parton system has the correct quantum numbers. We view such 
system as an "effective" heavy quark-antiquark pair that is color-singlet and this allows us to obtain approximate 
wavefunctions for the higher Fock components. To our knowledge exact solutions for the n > 3 Fock components of 
mesons do not exist. 

B. The light-cone wavefunction 

Simulations of high-p^ quarkonium dissociation in the QGP require knowledge of the light-cone wavefunctions of 
the state, which can be represented as follows [321 [33j: 

oo „ n ,2, , In \ / n 

|P+;J) = 4(P+;J)|0)= £ / II l 3 7 X ^ I*! «0 * E *i ~ 1 W E kj 

x I • • • al(x qt P+ + k qt , a qi ) ■ ■ ■ bl( X g j P+ + k^. )C% ) 4 k (x gk P+ + k 3fc ,a 9 J • • • ) . (23) 

Here, P+ = (P + ,P) are the large light-cone momentum and transverse momentum components of the hadron, and 
the momenta of the partons are given by {xiP + ,xtP +k,). In Eq. ()23[) a Qi is a set of additional relevant quantum 
numbers, such as hclicity and color. From the normalization of the meson state of fixed projection A of the total 
angular momentum J: 

(P+; J\P+'; J) = 2P+{2tt) 3 S(P + - P+')5 2 (P - P')S X \> , (24) 

an integral constraint on the norm of the light-cone wavefunctions tp(xi, kjj a{) can be derived. 

In the context of Eq. (|23p , the color-singlet contribution to quarkonium production can be understood as one that 
is already in the dominant, lowest order (n = 2), Fock component in the hard scattering. A color-octet contribution 
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requires the emission of at least one gluon for a color neutral hadron to be produced. In cither case the hadron state 
can be approximated as: 

\P + ) = J ( j^ a 2 . ^_ ^f *(?, k)a^ ci (xP + + k)b^((l - x)P + - k)|0) , (25) 

where for color-octet states a^ (6^) represent an "effective" heavy quark (anti-quark) in the 3 (3) state and [33[ 

*(x, k) = Norm x exp (- J^j^ j , ^ |d*d 2 k kM*,k)| 2 = 1 . (26) 

Note that in Eq. (|26[) we have used the fact that for quarkonia (1 — x)rriQ + xttlq = mg. If we introduce the notation 
Ak = ki — k2 — 2k, the transverse width A of the light-cone wavefunction ip(x, k) is determined from the condition: 

^L^jdxd'k Ak 2 |^,k)| 2 = 4(k 2 ) = \k> . (27) 

The factor 2/3 comes from the 2D projection of the mean squared transverse momentum k 2 from the instant wave- 
function form calculated in Tables U HU 

V. QUARKONIUM DYNAMICS AT HIGH TRANSVERSE MOMENTUM 

In this section we present the details of the dissociation model of quarkonium propagation through the QGP. The 
essence of the dissociation model for heavy mesons is that they have short formation times and can therefore form in 
the medium on a time scale iform- Interaction with the thermal medium can dissociate the mesons on a time scale 
tdiss- The final yield is given by a rate equation which takes into account the formation and dissociation processes. 
In the next section we first discuss the rate equation abstractly, using if or m and tdiss as parameters. In the later 
sections we will estimate iform and calculate foiss- For more details on the dissociation model and its application to 
the phenomenology of open heavy flavor, see |33J ] . 

A. The rate equations 

Let us denote by N^= ) d (pT, a) the number of perturbatively produced point-like QQ states at transverse momentum 
Pt- Up to an overall multiplicative Glauber scaling factor Tab, these pr number distributions are directly proportional 
to the cross sections discussed in Section ITT] The time scale for the hard QCD process is given by t ~ 1/rriT, where 
rriT = \/p\ + M 2 . For transverse momenta above a few GeV and M c5 > 3 GeV, M b i > 9 GeV, respectively, 
this production time is very short. Thus, we take this time to be the starting point for the evolution of the QQ 
state. The rate of formation of the corresponding hadronic state, with the appropriate quantum numbers {a}, is 
given by the inverse formation time l/if or m(PT)Qf). In the presence of a medium, the meson multiplicity, which we 
denote by N^ on (p T ,a), is reduced by collisional dissociation processes at a rate l/idiss(PT,aO- Finally, the number 
of dissociated QQ pairs with a net transverse momentum pt is N dl ^-(pr,a). At the transverse momenta that we 
consider, the probability for the heavy quark Q or antiquark Q will pick up a thermal partner and reform a quarkonium 
state is negligible. The heavy (anti)quark fragmentation contribution to quarkonia is also negligible. Heavy quarks 
fragment primarily into open heavy flavor mesons. 

The dynamics of such a system is governed by the following set of ordinary differential equations: 

m = ~7 u \^QQ ^'PT, a )i (2») 

at tf OTm (t;p T ,a) w 

QQ ;N^A d (t;p T ,a)-- rN%^(t; PT ,a), (29) 



dt tfoim(t;pT,a) QQ ' ' tdi SS .(t;p T ,a) 

N^ son (t- lP T,a), (30) 



dN*% s -(t;p T ,a) i 



dt t diss (t;p T ,a)' 
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subject to the constraint N^ d (t;p T , a)+N^ s ' 
by the initial conditions 



\t- lPT ,a)+N^g-(t ]PT ,a) 



QQ 



Nq% (j>t, «), and is uniquely determined 



N%% d (t = 0;p T ,a) 

N^ s ™(t = 0;p T ,a) 

N«% s (t = 0;p T ,a) 



N^A d (pQCD; PT ,a), 



"QQ 



0, 
0. 



(31) 
(32) 

(33) 



Note that in Eqs. ((30|) the evolution of the dissociated QQ pair into D- or B- mesons is not shown since it does not 
couple back to Eqs. (|25j) . (|2^1) . 

Realistic simulations include the velocity and {a} dependence of the formation rate of all quarkonium states and 
the time, position and {a} dependence of their dissociation rate. It is, however, useful to integrate the system of 
equations analytically for a simple test case. Our simplified test case assumes that the dissociation time is constant 
in the interval < t < L and if t > L. The solution for the QQ mesons as a function of time is: 



N^ son (Q<t<L;p T ,a) 



"QQ 



Armcson 

^QQ 



(t > L;p T ,a) 



N^ d (pQCD- PT ,a) 



N^ d (pQCD- PT ,a) 



= -L/tf orm (pT,a) 



tdi ss .(PT,a) 



tdiss. {PT, Ot) - t forT:1 (pT , a) 
tdiss.(PT,a) 



^diss. (pt, a) - if orm (p T , a) 

-t/tform (PT, a) \ 



( 



diBs.(Pr,a) 



"t/tform (PT, a) 



,(34) 



.(PT,a) _ e -L/t form (p 



T,a) \ 



(35) 



The interested reader can easily deduce the solutions for N^ d (t; pt , a) and N dl g-(t;pT,a) and verify that the 

solutions in Eqs. ([34| . (|35j) are finite for tf OTnl (pT,a) = idiss. (pt, «)■ Eqs. (l34|) . |35|) can be used to understand the 
qualitative features of the time dependence of quarkonium formation. 



B. Formation time of quarkonium states 



The approach to estimating the formation time of quarkonium states differs considerably from the approach used 
for open heavy flavor J32tl33l| or light particles 55j that come from the fragmentation of a hard parton. In the latter 
case the formation time is inversely proportional to the virtuality of the parton decay and is governed by longitudinal 
dynamics. For quarkonia, the QQ state is prepared instantly (~ 1/ \fpj- + M 2 ) in the hard collision and subsequently 
expands to the spatial extent determined by the size of the asymptotic wavefunction. In this case all spatial directions 
are important. The velocity of the heavy quarks in the meson and a typical upper limit of the meson formation time 
can be evaluated as follows: 



rest frame 



lh 



(36) 



where the typical momenta, k, and sizes, a, are given in Tables HI HTl In this paper we are interested in high transverse 
momentum mesons, in which case there is a boost in the direction of propagation and, consequently, time dilation 



C£(Pr,a)=7* 



rest frame 



(a) = 7-5- 
PQ 



1 = 



Vp 2 t + RP 
M 



(37) 



For example, for the relevant formation time determined by the expansion of the QQ state in a direction transverse 
to the direction of propagation the transverse size remains the same when boosted back to the laboratory frame but 
the velocity transforms by picking up a factor of 1/7. Note that in Eq. (I37p \p\ = pt since in this paper we work at 
mid-rapidity. The masses of the quarkonium states M are taken from [56j. 7 is the meson boost factor. Since the 
formation process is non-perturbative and can not be modeled accurately, the values of if orm obtained from Eq. [37] 
should be considered as an estimate of the upper bound. Therefore, in addition to calculating the final yields for 
^form = ifmax = op j we also calculate the yields for if or m = ifmin = 2§~ and the variation gives us an estimate of the 
uncertainty due to the uncertainty in the formation time. This does not change our results qualitatively. 
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Dissociation time of quarkonium states 



The propagation of a QQ state in matter is accompanied by collisional interactions mediated at the partonic level, 
as long as the momentum exchanges can resolve the partonic structure of the meson. Two effects are related to 
these interactions: a) a broadening in the distribution of quarkonium states relative to the original direction; b) 
a modification of the quarkonium wavcfunction. The former effect integrates out as long as we consider inclusive 
production. The latter effect leads to the dissociation of the meson state. 

Let us define: 



xfc 



dr 



M 2 (x(r),r) 

V x ( r )> r ) 



£: 



x(t) = x +/3(t-t ) 



(38) 



In Eq. 



(jl 2 is the typical squared transverse momentum transfer given by the Dcbye screening scale, [i = gT for a 



gluon-dominated plasma. X q is the mean scattering length of the quark and £ ~ few is an enhancement factor from 
the power law tail of the differential scattering cross section. Finally, Xo is the position of the propagating QQ and 
f3 is the velocity of the heavy meson. Note that \/3\ < 1. For a medium of uniform partem density and length L 
X/i 2 £ = /j, 2 (L/\ a )$. For the realistic expanding medium the partem density and temperature can be determined as 
described in |57l l58l ] and the integral in Eq. (|38[) can be performed numerically. 

With the results for the cumulative momentum transfer at hand, the medium-modified quarkonium wavcfunction 
can be evaluated analytically for the functional form specified in Eq. (|27p . The survival probability for the closed 
heavy flavor mesons is given by 



-Rmrv.(x^ 2 £) = 



1 



2(2tt) 3 

1 



tfkdx i>) ( Ak, x)i\)i ( Ak, x) 



2(2tt) s 



dx Norm itx{1 - x)A z e *<i—oa* 



2^(1 - x)AVx^ 2 C + x(l - x)A 2 



yjx{\ - x)A 2 + v /%^ 2 £ + x(l - x)A 2 



The dissociation rate is then then given by 



tdiss. (pt, a) 



dPm 



diss. 



dP a , 



dt 



dt 



(39) 



(40) 



The uncertainty in tdi ss arises from the uncertainty in the coupling between the heavy quarks and the medium 
(described by the strong coupling constant g) and the enhancement that arises from the power law tails of the Molicrc 
multiple scattering in the Gaussian approximation to transverse momentum diffusion |32j (described by £). In our 
calculations we use two sets g = 1.85, £ = 2 and g = 2, £ = 3 to understand the sensitivity of the results to these 
parameters. 

To get a sense of the formation and dissociation times involved in the quarkonium dynamics in the QGP, we give 
the values for a fixed transverse momentum pt = 10 GeV for the quarkonium states and consider their production in 
the center of the collision geometry. The specific nuclear collisions that we chose are 0-20% central Au+Au collisions 
at RHIC at \^Snn = 0.2 TcV. In the examples that follow g — 2 and £ = 2. Since the medium expands after the 
collision, we quote the formation and dissociation times at t = 1.5 fm/c. Results for charmonia and bottomonia are 
presented in Table III and Table IV, respectively. 



Charmonium state 



Formation time max [fm/c] 
Dissociation time [fm/c] 



J/iP 



X=0,l,2 



3.49 
1.73 



3.35 
1.74 



4.40 
1.61 



TABLE III: Upper limit on the formation time and dissociation time of quarkonium states of pr = 10 GeV and produced in 
the center of the nuclear overlap region of 0-20% central Au+Au collision at RHIC. 



Bottomonium state 


vtW 


T(l) 


»»(2) 


T(2) 


%(3) 


T(3) 


XM),1,2(1) 


Xmi,i.2(2) 


Xw,i,2(3) 


Formation time max [fm/c] 
Dissociation time [fm/c] 


1.44 
3.30 


1.44 
3.30 


2.84 
2.23 


2.85 
2.23 


4.24 
1.92 


4.17 
1.93 


2.36 
1.93 


3.45 
2.06 


6.23 
1.73 



TABLE IV: Upper limit on the formation time and dissociation time of bottomonium states of pt 
the center of the nuclear overlap region of 0-20% central Au+Au collision at RHIC. 



10 GeV and produced in 
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FIG. 3: Theoretical model predictions for J/ijj Raa in 0-20% central nucleus-nucleus collisions. The top left panel is for RHIC 
Au+Au collisions at \/£ = 0.2 TeV HH. The top right panel is for RHIC Cu+Cu collisions at ^5 = 0.2 TeV [||. The lower 
panels are for LHC Pb+Pb collisions at y/S = 2.76 TeV [6y,|6l|. The various curves represent the uncertainty due to variations 
in g and £ (taking the sets g — 1.85, £ = 2 and g = 2, £ = 3), and a factor of 2 variation in £f orm . These results exclude the 
Cronin effect. Dashed curves include the B — >• J/tp feed-down. 



VI. NUMERICAL RESULTS FOR THE NUCLEAR MODIFICATION FACTORS 



A. Raa for J/ip and T(15) 



In this section we exclude the Cronin effect but include initial-state cold nuclear matter energy loss and shadowing. 
Let us first consider the nuclear modification factor for J/tp mesons. From Fig. [3] we see that the measurement of Raa 
in RHIC Au+Au collisions [59| shows a suppression factor of about 0.8 at pr ~ 6 GeV. On the other hand, Cu+Cu 
collisions [49( show an enhancement factor of nearly 1.5. Both measurements have large error bars. Even including the 
uncertainty in our model parameters, we obtain a somewhat higher suppression (for pt ~ 6 GeV, Raa ~ 0.35 — 0.45 
for Au+Au and Raa ~ 0.45 — 0.65 for Cu+Cu) than the one currently observed at RHIC. In Fig. [31 our results for 
the prompt yields of J/tp mesons are marked by upper and lower yellow bands corresponding to the upper (if^m) 
and lower (ij-^m = C/^) limits of our formation time estimate respectively. The bands themselves correspond to 
our estimate of the uncertainly in the sets of parameters that determine the coupling of the heavy quarks with the 
in-medium partons [g = 1.85, £ = 2 (minimum considered coupling gives the upper limit of the yellow band) and 
g = 2, £ = 3 (maximum considered coupling for the lower limit of the yellow band)]. The pronounced effect of the 
variation of the formation time can be intuitively seen as follows. From Eg. 1301 we see that the dissociation mechanism 
is operative only when N^q 071 is substantial, i.e. after t( oml . Since the upper limit for formation time of quarkonia 
can be on the order of several fm/c, (see Tables Hill IIV[) . the density of the medium at if^rm ^ s reduced considerably 
due to Bjorken expansion, giving weaker dissociation and smaller suppression. This effect is more pronounced than 
the details of the coupling of heavy quarks to the in-medium partons. The RHIC experiments report suppression for 
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FIG. 4: (Color online) Expected J/ip suppression versus centrality (-/V par t) at \ // ~S = 2.76 TeV. The left panel compares the 
theoretical results to the ATLAS Rev measurement [61|]. The right panel compares the theoretical results to the CMS Raa 
measurement I60II. 



the inclusive J/ip yield. For direct comparisons, we also show the Raa for the inclusive yields in Fig. [3] with dashed 
lines. The color scheme for the various parameters is analogous to the direct production, though we do not color in 
yellow the band associated with the uncertainty in coupling to avoid cluttering. The B— meson yields for p+p and 
A+ A collisions were taken from [33[ . 

The high-pr suppression of J/ip mesons in Pb+Pb collisions at the LHC, reported by the ATLAS and CMS 
experiment, is substantially higher than the RHIC results. For example, for pt G (6.5, 30) GeV, CMS [601 shows 
Raa ~ 0.2 for the most central collisions. For a sharply falling spectrum we expect the suppression to be dominated 
by the low momenta in the px bin. (The mean px of the observed particles in the bin is roughly 9.3GcV, for 
p+p collisions.) On general grounds, one expects the suppression at LHC to be larger than at RHIC. The initial 
temperature at LHC is higher than at RHIC, and consequently the gluon density at any given time in the Bjorkcn 
expansion is also higher. This will give rise to more rapid dissociation at LHC both in equilibrium and non-equilibrium 
approaches. (The system size for Au and Pb are roughly the same.) 

For px = 9 GeV at the LHC, our calculations give Raa ~ 0.35 — 0.45, which is slightly smaller than at RHIC 
(at the same pr) but underestimates the observed suppression at LHC in the most central Pb+Pb collisions. The 
ATLAS experiment has presented the ratio of binary collision scaled central-to-peripheral J/ip yields Rcp [61j . More 
specifically, their baseline in given by the 40-80% peripheral Pb+Pb reactions. The ATLAS results are dominated by 
Pt > 6.5 GeV. (The cuts are such that 80% of the yields come from pt > 6.5GeV.) Comparison to the theoretical 
model calculation is shown in the left bottom panel of Fig. [31 The two data points are for 0-10% and 10-20% 
centrality. The CMS prompt J/tp Raa result [60| is shown in right bottom panel of Fig. [3j The two data points are 
again for 0-10% and 10-20% centrality and exhibit larger suppression than the theoretical model predictions. The 
ALICE experiment at LHC has measured the nuclear modification of J/tp production at forward rapidity [621 ] and 
low transverse momentum. To constrain theoretical models of quarkonium production in heavy ion collisions, it will 
be very helpful to extend these forward measurements to high transverse momentum. 

In our formalism, we have approximated the quarkonium wavefunction by the vacuum wavefunction, which is valid 
if the thermal effects on the quarkonium wavefunctions are small. The thermal wavefunctions (for example, those 
obtained by solving the Schrodinger equation in the thermal potentials @, LDtSI) will be wider in position space at 
higher temperature and therefore will dissociate more easily. Therefore a higher suppression at LHC could be the 
evidence for thermalization effects at the level of the quarkonium wavefunction at LHC. We leave a more detailed 
analysis of thermal effects on the wavefunctions for future work. Nevertheless, it is important to identify at what 
centrality the discrepancy between the present theoretical model predictions and the data appear. Fig. [4] we show the 
PT-averaged suppression, 



RAA{N palt ) (or R C p(N part )) 



J . dp T Raa {pt ', -/Vpart ) (or R C p {pt ; N paxt ) ) 



da 
dydpq 



L ra d PT 



(41) 



dydpq 



of J/ip mesons versus centrality. We present a comparison to the ATLAS central-to-peripheral data KJlf in the left 
panel. The deviation between data and theory is only for N pajrt > 300. A comparison to the CMS data [6(| is shown 
in the right panel. In this case the deviation between data and theory is for N part > 200. 
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FIG. 5: (Color online) Theoretical model predictions for T Raa in 0-20% central nucleus-nucleus collisions. The top left panel 
shows Raa results for T(15 r ) + T(2S) + T(35) in Au+Au at V5 = 0.2 TeV. The bottom left and right panels shows Raa 
results for T(2S) + TJ3S) and T(IS'), respectively, in Pb+Pb at y/S = 2.76 TeV. The top right panel data from CMS is for 
— 100% centrality [60] compared to the theoretical prediction at the minimum bias iVpart ~ 110. 



The CMS experiment at the LHC has also measured Raa for T(nS) states in Pb+Pb collisions at \/S — 2.76 TeV 
per nucleon pair. The result [63[ is presented for decay muons satisfying the transverse momentum cut Pt(i^) > 4 GeV 
and rapidity |?y| < 2.4. 



T(25 + 3S). 



T(1S) 
T(2S + 3S) 
T(15) 



0.76I^;^±0.12, 



(42) 



PbPb = 0.24i°; x J ± 0.02 , 



giving 



R AA (T(2S + 3S)) 
R AA (T(1S)) 



0.32to;is±0.03 



(43) 



We cannot calculate the equivalent ratio of RaaS because our formalism for the production and propagation of 
Ts is not applicable to pt(Y) < 6 GeV. In our approach the meson should be boosted relative to the medium. 
Furthermore, for static or slowly moving mesons one may expect thermal effects on the quarkonium wavefunction. 
These are precisely the Ts that determine the total yield ratios in Eqs. (|42"T) . (|43|) . 

For pr > 7 GeV at the LHC, our calculations show that j? <r(is\) ~ ^ (Fig- E]). Physically, this is because 
the higher dissociation rate for the excited mesons is compensated by their larger formation time, which means 
that dissociation becomes dominant after the medium has diluted. If in future experiments a p^-differcntial yield of 
T(2S' + 35) shows a much larger suppression compared to T(15') for p? > 7 GeV, that may also support the possibility 
for thermal effects on the quarkonium wavefunction at high px- This is because the thermal wavefunctions of higher 
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excited states, being broader in position space, are more easily affected by the QGP, thereby naturally giving even 
higher dissociation rates. Fig.[S]also shows theoretical predictions for the suppression of the various T states in central 
Au+Au and Pb+Pb collisions at RHIC and the LHC, respectively. In the top left panel present a comparison to the 
minimum bias pr-differential T(15 l ) CMS nuclear modification measurement. In this case, the theoretical calculation 
is performed for the average number of participants for minimum bias collisions. 

The overall suppression of the J/ip and the T yields in A+A collisions is a combination of CNM and QGP effects. 
In this section we ignored transverse momentum broadening effects. Therefore, the suppression is largely due to cold 
nuclear matter energy loss and QGP dissociation (the effect of shadowing is small). In Section TVI B[ we will show 
also results for p+A collisions where QGP effects are absent. (Sec Fig. [7] for J/ip and Fig. [8] for T.) For px ~ 6 GcV, 
RpA ~ 0-8 for both J/ip and T. Noting that in A+A collisions the CNM effects are amplified relative to p+A, we 
conclude that a significant part of the suppression in quarkonium yield in our calculation comes from cold nuclear 
matter energy loss. The situation is more complicated when we include transverse momentum broadening as we 
discuss next. 

B. Transverse momentum broadening effects 

As expected, transverse momentum broadening effects enhance the production of J/ip mesons for p? ~ 5 — 10 GeV 
in heavy ion collisions (Fig. [5]). The Cronin effect is not very important for pr > 10 GeV at these center-of-mass 
energies. The main features are as follows: the Cronin effect can substantially alter the Raa of J/ip mesons with 
Pt ~ 7 GeV at RHIC. For example, for Au+Au collisions Raa increases from ~ 0.35 — 0.45 to ~ 0.8 — 1.2. For 
Cu+Cu collisions Raa increases from ~ 0.45 — 0.65 to ~ 0.7 — 1.0. The reason behind this enhancement is that the 
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FIG. 6: (Color online) Theoretically calculated Raa for J/ip, including the Cronin effect. The top left panel is for RHIC 
Au+Au collisions at vS = 0.2 TeV. The top right panel is for RHIC Cu+Cu collisions at y/S = 0.2 TeV. The lower two panels 
are for LHC Pb+Pb collisions at \/S - 2.76 TeV as in Fig. El 
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mean scattering lengths of the initial-state gluons (which dominate quarkonium production) are considerably smaller 
than the mean scattering lengths for quarks. The Cronin effect at the LHC is smaller than the Cronin effect at 
RHIC because of the harder quarkonium spectrum. Including broadening only increases Raa at pt ~ 7 GcV from 
~ 0.3 — 0.4 to ~ 0.6 — 0.8. The theoretically calculated Raa for J/tp may be consistent with the RHIC Cu+Cu results 
but is only marginally compatible with the RHIC Au+Au results. Also, it is clearly incompatible with the ATLAS 
central-to-peripheral Pb+Pb measurement at the LHC and the CMS Raa suppression. The differences in the degree 
of pt — 5 — 10 GeV suppression of quarkonium production between RHIC and the LHC in Fig. [5] suggest that better 
understanding of the Cronin effect (if any) is necessary for consistent J/tp phenomenology in heavy ion collisions. 

More specifically, experimental quarkonium measurements in p+A or d+A collisions, where effects from the QGP 
are absent, are important to constrain the cold nuclear matter effect. Our theoretical predictions for the high-p^ J/tp 
nuclear modification factor R p a in such collisions are given in Fig. [7] The upper (solid red) curve includes the Cronin 
effect. The lower (solid blue) curve includes only power corrections and CNM energy loss. We finally point out that 
including the contribution from B decays reduces slightly the Cronin enhancement for px ~ 5 — 10 GeV. 

The combined results for Raa (0-20% central) and R p a (minimum-bias) for bottomonia are given in Fig. |8] The 
effect of transverse momentum broadening is much smaller for bottomonia when compared to the one for charmonia. 
This can be intuitively understood as follows. The mechanism for Cronin enhancement in this calculation is that 
initial state scattering increases the typical transverse momentum carried by the incident partons by a few GeV. For 
quarkonia, there is an additional scale mjy. For bottomonia the mass scale is considerably larger than the transverse 
momentum broadening scale and few additional GeV do not increase the yield significantly. Preliminary T suppression 
measurements are now available at RHIC |64[ . More differential pr measurements will shed light on the similarities 
and differences in the CNM effects at RHIC and at the LHC. 
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FIG. 7: (Color online) Theoretical predictions for J/tp R p a in minimum-bias collisions with and without the Cronin effect. In 
this plot we only show the prompt yields. The top left panel is for RHIC d+Au collisions at VS = 0.2 TeV. The top right 
panel is for RHIC d+Cu collisions at VS = 0.2 TeV. The lower panel is for LHC p+Pb collisions at VS = 2.76 TeV. 
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FIG. 8: (Color online) Theoretical model predictions for minimum-bias R p a(T) and central 0-20% Raa{Y). The left panel 
is for RHIC p+Au and Au+Au collisions at y/S = 0.2 TeV. Data is from STAR [HI The right panel is for LHC p+Pb and 
Pb+Pb collisions at VS = 2.76 TeV. Data is from CMS H. 



VII. CONCLUSIONS 



In summary, we carried out a detailed study of high transverse momentum quarkonium production and modification 
in heavy ion reactions at RHIC and at the LHC. We used a Heavy Quark Effective Theory approach to calculate the 
baseline quarkonium cross sections. We found that for J/-0 mesons the theoretically computed spectrum is slightly 
harder than the one observed in the experiment. For all T states (IS*, 25*, 3S) the agreement is within a factor of two 
when we consider both the TeVatron and the LHC data. In reactions with heavy nuclei, we presented theoretical 
model predictions for the nuclear modification of quarkonium yields at high px in minimum bias p(d)+A and 0-20% 
central A+A collisions. We focused on the consistent inclusion of both cold (CNM) and hot (QGP) nuclear matter 
effects in different colliding systems at different center-of-mass energies. We compared our results to published and 
preliminary experimental data, where applicable. 

In calculating the spectra of quarkonia in heavy ion reactions, we included nuclear shadowing (here implemented 
as coherent power corrections) and initial-state energy loss. We also provided, to the best of our knowledge, the 
first implementation of initial-state transverse momentum broadening to study phenomenologically a possible Cronin- 
like enhancement for quarkonia. All these effects have been well studied for light partons and have been recently 
incorporated in open heavy flavor production. In this paper we extended them to J/ip and T mesons. Since the QQ 
pair is created in the short-distance hard scattering process and evolves quickly into a component of the quarkonium 
wavefunction, the effects of propagation through the QGP were included through a collisional meson dissociation 
model, which was successful in describing the attenuation of open heavy flavor, B — >• £ + X, D — > £ + X, at RHIC. 
In this paper we restricted our results to high transverse momentum and explored the consequences of assuming that 
the initial wavefunctions of the quarkonia are well approximated by vacuum wavefunctions in the short period before 
the dissociation. 

We found that ignoring the Cronin effect leads to a small overestimate of the suppression of J/tp mesons in the px 
region between 5 GeV and 10 GeV in central Cu+Cu and Au+Au collisions at V^ = 0.2 TeV at RHIC. Including 
initial-state transverse momentum broadening leads to a somewhat better agreement between theory and the current 
experimental data only for the Cu+Cu reactions. A smaller Cronin enhancement will work better. We demonstrated 
that CNM effects can be easily constrained in d+A reactions at RHIC. For example, the d+Au calculation that 
includes power corrections and cold nuclear matter energy loss predicts 20% suppression of the J/tp cross section. 
Including transverse momentum broadening may lead to as much as 50% enhancement in the region of the Cronin 
peak. We also found that the Cronin-like modification of the T spectrum is much smaller. Current data on high-pr 
quarkonium production at RHIC does not indicate the presence of thermal effects at the level of the quarkonium 
wavefunction within our theoretical framework. 

The conclusions from our theoretical model comparison to the vS = 2.76 TeV LHC data are not as clear cut. Our 
calculations underestimated the suppression for J/ip production reported by the ATLAS and CMS experiments in 
the most central Pb+Pb collisions. On the other hand, it works quite well in mid-central and peripheral reactions. 
We found that the Cronin enhancement at the LHC was smaller than the one at RHIC due to the harder spectra. 
However, any Cronin enhancement appears incompatible with the experimental results. For T mesons, pr-diffcrcntial 
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data in A+A collisions is scarce. CMS measurements of minimum bias T(15) indicate that the low pt suppression 
may decrease or disappear at high p^. At the same time, at low p^, where our calculation is not applicable, CMS 
reported a large relative suppression of T(2S + 35) to T(IS'). If this measurement is extended to high pr with 
similar results, it will clearly be incompatible with our model predictions with quarkonium wavefunctions unaffected 
by thermal effects. Together with a more refined p^-differential measurement of J/ip suppression at yS = 2.76 TeV, 
this will be a strong indication that thermal QGP effects may persist for high transverse momentum quarkonia at the 
LHC in central Pb+Pb reactions. We plan to address this possibility in a separate publication. 
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Appendix A: Baseline for LHC at VS = 2.76 GeV 

Here we show (Fig. [5J the cocktail of contributions that give the baseline production yield for J/ip at LHC at 
y/S = 2.76 GeV in the HQET formalism. 

For the direct production, as described in Section HH there are three contributions: one color-singlet contribution 
([ 3 Si]i) three color-octet contributions ([ 1 5o]s), (P^ils); and ([ 3 Po]s)j with matrix elements given in Section [Til 

For the prompt yield, we add the feed-down contributions from the XcJ- The contributions of the XcJ are also 
shown in Fig. [9] To obtain the feed-down contributions, we multiply the corresponding yields by the corresponding 
branching fractions, 

BR(xco(l) ->- J/4') = 0.0116 , 

BR(xd(l) ^ JM = 0.344, (Al) 

BR(xc2(l) -► J/VO = 0.21 . 

This gives the p+p baseline contribution for the J/i/j- We have included experimental data on prompt J/tjj production 
at \/S = 2.76TeV [g3|. 

Each of the species will undergo a different modification due to CNM and QGP effects. To obtain the p+A and 
A+A yields, we calculate the modified yields for each species and combine them again to obtain the net modification. 
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